Strong magnetohydrodynamic turbulence with cross helicity 
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Magnetohydrodynamics (MHD) provides the simplest description of magnetic plasma turbulence 
in a variety of astrophysical and laboratory systems. MHD turbulence with nonzero cross helicity 
is often called imbalanced, as it implies that the energies of Alfven fluctuations propagating parallel 
and anti-parallel the background field are not equal. Recent analytical and numerical studies have 
revealed that at every scale, MHD turbulence consists of regions of positive and negative cross 
helicity, indicating that such turbulence is inherently locally imbalanced. In this paper, results 
from high resolution numerical simulations of steady-state incompressible MHD turbulence, with 
and without cross helicity are presented. It is argued that the inertial range scaling of the energy 
spectra (E^) of fluctuations moving in opposite directions is independent of the amount of cross- 
helicity. When cross helicity is nonzero, and E~ maintain the same scaling, but have differing 
amplitudes depending on the amount of cross-helicity. 

PACS numbers: 



I. INTRODUCTION 

Magnetohydrodynamic (MHD) equations govern the 
dynamics of a plasma when the spatial scales of inter- 
est greatly exceed any intrinsic plasma scale. The MHD 
model is relevant in understanding waves, instabilities, 
and turbulence in a variety of astrophysical systems. 
For instance, magnetohydrodynamic (MHD) turbulence 
plays a key role in the theoretical modeling of the spectra 
of velocity and magnetic fluctuations in the Solar Wind 
[e.g., as well as electron density fluctuations in the 
interstellar medium [e.g.,0,0]- 

Theoretical modeling of the spectra in strong MHD 
turbulence has developed on different fronts: phe- 
nomenology, statistical closures and more recently, high- 
resolution numerical simulations. The first phenomeno- 
logical model, a la Kolmogorov @, , of the energy spec- 
trum was developed independently by Iroshnikov 10] and 
Kraichnan [ll|. In these works they realized that the 
magnetic field of the large scale eddies acts in much the 
same way as a guide field that supports smaller scale 
Alfvenic fiuctuations, and the turbulent energy transfer 
takes place as the result of many cumulative collisions 
between counter-propagating Alfven wave packets trav- 
eling along the local magnetic field. One deficiency of 
this phenomenology is that it is based on the assumption 
of an isotropic spectral transfer, in contradiction with 
more recent results that reveal the anisotropic charac- 
ter of MHD turbulence [e.g., [l^, To account for 
anisotropy in the strong turbulence regime, Goldreich 
and Sridhar [ijj, hereafter referred as GS, proposed a 
new phenomenology in which they introduced the so- 
called critical balance condition, stating that the turbu- 
lence is strong when the time scale for wave propaga- 
tion along the magnetic field matches the nonlinear in- 
teraction time. In the past decade, constantly increasing 
massively parallel computer clusters have become large 
enough to allow one to simulate the inertial range spectra 



of strong MHD turbulence, leading to an extensive vol- 
ume of work addressing universal scaling laws in MHD 
turbulence [e.g., Il5l425l |. These simulations have moti- 
vated new phenomenological models and have resulted 
in a renewed interest in the fundamentals of MHD tur- 
bulence. 

One aspect of MHD turbulence, which has recently 
drawn significant attention is the role that cross-helicity 
plays in the turbulent cascade [25l - l30j . Both total en- 
ergy E = ^ J{v'^ + b'^)cPx and total cross-hehcity H'^ = 
J (v-h)(Px are conserved by nonlinear interactions; where 
V and b are velocity and magnetic fields. In MHD turbu- 
lence both energy and cross helicity cascade from large 
to small scales as a result of the nonlinear interaction 
of counter-propagating Alfven fiuctuations. If we denote 
= I /(v ± b)2 (Px, which are the energies carried 
by Alfven fiuctuations propagating in opposite directions 
(the details are given below) , the total energy and cross- 
helicity of the system are given hy E = E~^ + E~ and 
He = E'^ — E^ , respectively. Therefore, when the latter 
is nonzero, the energies of waves propagating along and 
against the guide field are not equal, and in that sense 
the turbulence is called imbalanced. 

Independent numerical simulations by different groups 
have demonstrated that strong MHD turbulence is al- 
ways locally imbalanced, in the sense that in the steady 
state, the turbulence develops correlated regions of posi- 
tive and negative cross-helicity, irrespective of the overall 
cross-helicity of the system. Imbalance turbulence is also 
present in the solar wind, where velocity and magnetic 
fluctuations show high correlations of a preferred sign, 
that is, the normalized cross-helicity ac = (v • b) /E = 
Hc/E is close to unity. The preferred positive sign of 
CTc indicates that there is more energy in Alfven waves 
propagating outwards from the Sun than propagating in- 
wards 3]. 

Recently, several phenomenological models have ad- 
dressed steady-state strong imbalanced MHD turbu- 
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lence [23|, I25l - l27l |29f, some with support from numeri- 
cal simulations and observations. However, these models 
have lead to conflicting predictions. For instance, the 
theories by Lithwick et al [26] and Beresnyak and Lazar- 
ian ^23] conclude that in imbalanced regions the Elsasser 
spectra have the same scalings E~^{k±) oc E^{k±_) oc 
^^5/3^ The theory by Chandran '2?] proposes that the 
spectra of E'^{k^) and E~{kj_) are different depending 
of the degree of imbalance. Finally, the analysis by Perez 
& Boldyrev [2^ and Podesta & Bhattacharjee [29| finds 
that the spectra of E~^{k±) and E^{k±) have different 
amplitudes but the same scalings E^{k±) cx E^(k±) cx 
I -3/2 
^ . 

In this work we present results from high resolution di- 
rect numerical simulations of strong and steadily driven 
MHD turbulence, aimed to elucidate the role of cross- 
helicity and resolve the controversies among different the- 
ories. Hereafter, strong turbulence is defined in the sense 
of Golreich and Sridhar [TJ]. This paper is organized 
as follows. Section [IT] briefly describes the MHD equa- 
tions and the relevance of using the Reduced MHD model 
for strong turbulence simulations. Section Hill describes 
in detail the proposed numerical strategy and discusses 
important numerical aspects that need to be considered 
when simulating strong imbalanced MHD turbulence. In 
section IIVI the most relevant results from an extensive 
number of numerical simulations are presented. In sec- 
tion |V] a phenomenological model for the simulation re- 
sults is proposed and section IVll presents the conclusions. 



II. MODEL EQUATIONS 

In the presence of a guide field Bq, the incompressible 
MHD equations describing the evolution of magnetic and 
velocity fluctuations, b(x,t) and v(x, t), can be written 
in terms of the so-called Elsasser variables [SfJ, = 
v±b: 



teractions take place when the fields overlap or "collide" 
with each other. It is convenient to describe the resulting 
energy spectra in terms of the so-called Elsasser energies, 
defined as E^ = J |z±|2d3x/4. 

In the limit of small amplitude of fluctuations, the in- 
compressible MHD system ([T|) describes non-interacting 
linear Alfven waves with dispersion relation a;^(k) = 
±fc||Wyi. The incompressibility condition requires these 
waves to be transverse, and they are typically decom- 
posed into two polarizations, the shear- Alfven wave (zg) 
and the pseudo- Alfven waves (zp) given in Fourier space 
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where ey is the unit vector in the direction of the guide 
field. Strong MHD turbulence is dominated by fluctua- 



tions with k^ ^ Goldreich and Sridhar argued 
that since for large k± the polarization of the pseudo- 
Alfven fluctuations is almost parallel to the guide field, 
such fluctuations are coupled only to field-parallel gra- 
dients, which are small since fc|| <C kj_. Therefore, the 
pseudo-Alfven modes do not play a dynamically essen- 
tial role in the turbulent cascade. We can remove the 
pseudo- Alfven modes by setting zj|^ = in equations ^ 
to obtain 

atz±T(VA-V)i± + (iT-V)i± = -ViP+iv2z±, (3) 

where R denotes the Reynolds numbers (discussed be- 
low). In this system, the fluctuating fields have only two 
vector components, = {zf, 0} (where we chose the 
z axis along the guide field Bq) but depend on all three 
spatial coordinates. It can be demonstrated that system 
([3]) is equivalent to the Reduced MHD model, originally 
developed for tokamak plasmas by [l^l and [33| . 



III. NUMERICAL APPROACH 
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where va — Bo/\/47rp is the Alfven velocity, p is the 
fiuid density, P is the total pressure that is determined 
from the incompressibility condition, V • z^ = 0, are 
large-scale forcing, and ly is the viscosity which acts as 
an energy sink at small scales. In contrast with a uni- 
form flow in hydrodynamic, which can be removed by 
a Galilean transform, the guide magnetic field, whether 
imposed externally or created by large-scale fiuctuations, 
is essential in MHD turbulence. It cannot be removed by 
Galilean transform and it mediates nonlinear interactions 
at all smaller scales. The linear terms, (v^ • V)z^, de- 
scribe advection of Alfven wave packets along the guide 
field, while the nonlinear interaction terms, (z^ • V) z^, 
are responsible for energy redistribution over scales. It 
can be observed from equations ([T|) that nonlinear inter- 
actions can only occur between z+ and z~, and such in- 



The reduced MHD model ^ describing the nonlinear 
interactions of Shear- Alfven modes is a valuable tool in 
theoretical and numerical studies of incompressible MHD 
turbulence. Based on this model, we now discuss the 
conditions that should be satisfied in order to correctly 
simulate strong MHD turbulence, both balanced and im- 
balanced. 



A. Computational domain 

In the presence of a strong guide field, MHD turbu- 
lence is inherently anisotropic. It is important to point 
out that such anisotropy will develop deep in the inertial 
range even if it is not present at the outer scale where the 
turbulence is driven. This can be understood from the 
Golreich and Sridhar [l3| picture of strong turbulence. 
As the turbulence cascades to smaller scales, eddies be- 
come more and more shrunk in the field-perpendicular 
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direction until fluctuations satisfy the critical balance 
condition fc||i?o ^ k±b\, where ~ l/l and k± ~ 1/A 
are the field-parallel and field-perpendicular wave vec- 
tors associated with an anisotropic eddy of parallel size 
I and perpendicular size A, respectively, and bx denoted 
rms magnetic fluctuations at scale A. It is assumed that 
bx ^ vx. 

In a traditional cubic simulation box, when an 
anisotropic wave packet fits the field-parallel direction 
(z— dimension of the box) its field-perpendicular dimen- 
sions are much smaller than the x and y dimensions of 
the box. This is obviously not an optimal situation, 
since the field-perpendicular resolution required to re- 
solve such a wave packet should be much higher than 
the field-parallel resolution. This decreases the effec- 
tive field-perpendicular resolution of the simulations. For 
example, in the case a cubic 512'^ box with the guide 
field Bq — 5brms, one will have an equivalent field- 
perpendicular resolution of only 512/5 ^ 100. 

For our simulations, we define the nonlinear interaction 
strength parameter 

X- (fc±6A)/(fc||Bo), (4) 

so that the critical balance condition then implies x ^ 1. 

Simulations of steadily-driven (dissipative) incom- 
pressible turbulence are generally based on the Fourier 
pseudo-spectral method. As with any spectral method, 
the solution to the differential equation is approximated 
by a truncated Fourier expansion, and the partial differ- 
ential equation is converted to a set of ordinary differen- 
tial equations in time for the N Fourier coefficients. Non- 
linear terms in this representation become convolutions 
whose direct computation requires 0{N^) operations. In 
pseudo-spectral methods, the convolutions are computed 
in real (or configuration) space by means of a Fast Fourier 
Transform (FFT) that only requires 0{NlogN) opera- 
tions. For the numerical simulations presented in this 
work, a third order semi-implicit Runge-Kutta/Crank- 
Nicholson method was used for the time integration. 

In order to allow for the incrtial interval to develop, 
turbulence is driven at the lowest resolvable wave num- 
bers, and the energy dissipates at large wave numbers 
determined by the Reynolds numbers. For simulations 
on a cubic periodic box of size L, the smallest wave- 
numbers along the field-parallel and field-perpendicular 
directions coincide, i.e., fcj^ = fc|| = 2t:/L. Therefore, 
driving at the low k\\ , k± results in an isotropic forcing 
and the nonlinear strength parameter at the forcing scale 
becomes xo = k±bx/k^\Bo ^ bx/Bo ^ 1, which means 
that at least at the large scales, nonlinear interactions are 
weak. This would not be harmful if we had the resources 
to achieve arbitrarily high resolution, as the turbulence 
would proceed weakly until x '^^^ Ij then it would 
become strong. However, simulations generally produce 
a rather limited inertial range, so that the parameter x 
can hardly reach unit y in such a setup. 

As pointed out by [la|, and applied in recent simula- 
tions pol . [21I [23 | , an effective way to avoid this is to use 



an anisotropic domain such that at the forcing scale the 
parameter x is already of order unity, that is, the ex- 
cited large-scale modes are already anisotropic and sat- 
isfy the critical balance condition. To achieve this we 
choose an elongated box x Ly, so that the lowest 
field-perpendicular and field-parallel wave-numbers are 
k± — 2tt/L± and fc|| = 27r/L||, respectively. In this case, 
forcing at the lowest k±,k\\ leads to xo = {L\\bx)/{L±Bo), 
which is of order unity provided that 

^ bx/Bo. (5) 

In this way, the turbulence is excited in a strong regime 
and the cascade proceeds down to smaller scales preserv- 
ing the critical balance condition. 



B. Numerical Resolution 

At first sight, it appears that elongating the box along 
the z direction to match the elongation of the eddies 
should not change the number of grid points required 
in this direction compared to the number of points in the 
X direction. Fortunately, the number of points in the z 
direction can be reduced. This follows from the fact that 
the turbulent spectrum declines quite slowly, as a power- 
law, in the k± direction, while it drops sharply in the fcy 
direction for fcy > fc ? . w here a is a some positive power 
not exceeding 1 [11,111, [H, HI [11] . This qualitatively 
different spectral behavior in fcy and k± directions allows 
one to reduce the numerical resolution by a factor of 2 to 
4 in the parallel direction, see Table HI We checked that 
the restoration of the full resolution in the z direction 
does not change the results, while significantly increases 
the computing costs. 



C. Periodic Boundary Conditions 

The spectral method assumes periodic boundary con- 
ditions in all spatial directions. The periodic bound- 
ary conditions in the direction of the guide field (the 
2— direction) may raise the question of whether the mag- 
netic field lines are periodic in such numerical simu- 
lations. If they were periodic then the Alfven modes 
counter-propagating along a given magnetic field line 
would repeatedly interact only among themselves, and, 
therefore, they might not become sufficiently decorre- 
lated between the consequent interactions. To answer 
this question we note that the numerical setup ensures 
periodicity of the fluctuations b, but not the magnetic 
field lines. Magnetic field lines are integral lines of a 
magnetic field, therefore, they are generally not peri- 
odic in the z— direction, since b(ky — 0) is generally 
nonzero. Consequently, an eddy traveling along a mag- 
netic field line interacts with many independent counter- 
propagating eddies. Reducing the parallel box size Ly 
may however lead to artificial overlap of an elongated 
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eddy with itself. This may explain why reducing the 
parallel box size Ly below ([5]) somewhat spoils the spec- 
trum at low wave numbers, as seen in forced simulations 
of Miiller & Grappin [loj. Increasing L|| beyond ([5]) does 
not change the results but increases the computational 
cost [Mason & Cattaneo, private communication 2006]. 



D. The strength of the guide field 

In the incrtial interval of turbulence, the guide field 
should be strong compared to the turbulent fiuctuations, 
b\ <^ Bq. It is important to establish how small the 
fiuctuations should be in order to exhibit the universal 
turbulent spectrum. The question is especially relevant 
for numerical simulations, as the large guide field implies 
small Alfvcnic time and, therefore, large computing cost. 
This problem was numerically addressed in [2(|, where 
it was found that the transition to the universal regime 
occurs approximately at bx/Bo ^ 1/3. For b\/Bo > 1/3 
the energy spectrum is closer to —5/3, possibly indicating 
that the magnetic field does not qualitatively change the 
Kolmogorov dynamics for the scales attainable in the nu- 
merical simulations. For bx/Bg < 1/3, the guide field sig- 
nificantly affects the dynamics, and the spectrum changes 
to —3/2. The "sweet spot" often used in numerical sim- 
ulations is brms/Bo ~ 1/5; it has been checked that the 
smaller ratio, bx/Bo ~ 1/10, does not lead to noticeably 
different results (see for example p^l - l2ll |). 

It is also important to note that in numerical simu- 
lations, where inertial intervals have quite limited ex- 
tent, the condition brms/Bo < 1/3 should be satis- 
fied at the outer scale. Indeed, since bx decreases with 
the scale quite slowly, say b\ cx A^/'^, if the condition 
brms/Bo < 1/3 is not satisfied at the outer scale, it can 
hardly be satisfied in the inertial interval. 



E. Reynolds numbers 

Probably the most significant limitation on numeri- 
cal simulations is imposed by the consideration of im- 
balanced turbulence. Indeed, in the imbalanced case, 
7 = / z~ > 1, the formally constructed Reynolds 
numbers Re^ = Xz^/v corresponding to z+ and z 
fields at some scale A would be essentially different (here 
z^ = One can argue that the effective Re number 

is then the smaller of the two, see below. Therefore, the 
resolution requirements increase with the amount of im- 
balance in order to produce large inertial ranges. Assume 
that the number of grid points in a field-perpendicular 
direction scales with the Reynolds number as iV ~ Re^, 
where /3 = 2/3 or 3/4 depending on the spectral slope 
(3/2 or 5/3). Then increasing the imbalance 7 by, say, 
a factor of 3 will require increasing the resolution by ap- 
proximately a factor of 2. Noting that resolution of at 
least 1024 is required to simulate the imbalance 7^2 
(see below), we conclude that significantly stronger im- 



Field-parallel energy spectrum at mode 




FIG. 1: Field-parallel energy spectrum at the dominant mode 
k± — 3 for different forcing correlation times. The dotted 
line represents the normalized spectrum of the forcing. We 
observe that as the correlation time increases, large fcy modes 
get suppressed despite the fact that force amplitude is the 
same for all modes. 

balance is not achievable with present day computing 
power. 



F. Random Forcing 

Another question that arises in the strong turbulence 
regime concerns the type of forcing. In real systems, 
large-scale turbulence can be driven by instabilities ex- 
citing both velocity and magnetic harmonics, by exter- 
nal antennae exciting magnetic fields, etc.; this should 
not matter for the spectrum of small-scale fluctuations. 
Similarly, in numerical simulations one can force at large 
scales either velocity or magnetic field, or both of them, 
this does not change the result ^21(1 . The goal is therefore 
not to mimic any real-life driving but rather to optimize 
the transition to the inertial interval. 

In this section we discuss the important aspects to be 
considered when choosing a particular forcing. We as- 
sume a random force f that has no component along z, it 
is solenoidal in the x — y plane and its Fourier coefficients 
outside the range 

1 < fc_L < 2, {2tt/L\\) < k\\ < {2TT/L\\)n^ (6) 

are zero, where determines the width of the force spec- 
trum in and Lj_ = 2tt. The Fourier coefficients inside 
that range are Gaussian random numbers with ampli- 
tude chosen so that the resulting rms velocity fluctua- 
tions are of order unity. The individual random values 
are refreshed independently at time intervals r. The pa- 
rameter riz controls the degree to which the critical bal- 
ance condition is satisfled at the forcing scale. Note that 
we do not drive the fc|| = mode but allow it to be gen- 
erated by nonlinear interactions. 

In contrast with incompressible hydrodynamic sys- 
tem, an incompressible MHD system can support Alfven 
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waves. When the system is driven by a time depen- 
dent forcing, the most effectively driven modes are those 
resonating with the frequencies present in the forcing. 
Therefore, the spatial spectra of the large-scale velocity 
and magnetic fluctuations are not generally the same as 
the spatial spectrum of the force. Rather, they essen- 
tially depend on the both spatial and temporal spectra 
of the random forcing. Ultimately, it is the spectrum 
of the large-scale velocity and magnetic fluctuations, not 
the driving force, that should be controlled in numerical 
simulations. To illustrate this we performed a series of 
simulations in which we drive the modes given by equa- 
tion ([6]) with Hz — 8. We would expect this to excite all 
modes from fco|| — {2tt /Lz) to 8fco||. Figure [T] shows the 
field-parallel energy spectrum of fluctuations at /cj. = 3, 
for differing values of the force correlation time, show- 
ing that for short time correlations, the fluid response 
broadens. 

This adds additional complexity as the way in which 
the turbulence is forced may affect the outcome of 
the simulation's results. In MHD turbulence there are 
regimes where the spectrum is expected to obey univer- 
sal power laws, as well as transition regions where the 
turbulence changes character, for instance from weak to 
strong turbulence. Therefore, when simulating MHD tur- 
bulence, meaningless results could easily be obtained if 
the forcing is not chosen carefully. Not optimized forcing 
makes the transition from the forcing scale to the iner- 
tial range unnecessarily longer, as the turbulent fluctua- 
tions have to develop the anisotropy implied by critical 
balance. In the following we perform simulations with 
a short-time-correlated random forcing that drives the 
turbulence close to the critical balance condition at the 
large scales. A short-time correlated Gaussian random 
force has another important advantage, that it allows 
one to control the rates at which the energies of z'^ and 
z~ modes are injected. 



G. Inertial range and bottleneck 

In spite of the significant growth in massively paral- 
lel supercomputers that has occurred over the last years, 
simulations of MHD turbulence still result in a rather lim- 
ited inertial range. This generally raises the question of 
as to whether the simulations fully capture all scales from 
forcing to dissipation and whether the scaling exponents 
inferred from the energy spectrum are accurate enough 
to confirm or rule out the existing models of turbulence. 
In order to extend the inertial range, some groups have 
performed simulations with hyper-viscosity of different 
orders. However, using artificial viscosity might enhance 
possible bottlenecks (or bumps) that arise at high k due 
to an abrupt viscous suppression of the turbulent cascade 
in the dissipation range, e.g., [H, [sl]. This bottleneck 
can significantly affect measured scaling laws in simula- 
tions with inertial ranges of limited extent. We avoid 
the use of hyper-viscosity and perform direct numerical 



simulations of the RMHD equations ^ . 

In our results we identify the inertial range by per- 
forming a set of simulations from low to high Reynolds 
numbers (with increasing resolution). It is verified that 
the inertial range becomes larger as the Reynolds number 
increases. At larger k the inertial interval is followed by 
the dissipation range where the power-law scaling does 
not apply. 

Spectral power laws in the simulations are detected 
by compensating the energy spectra with corresponding 
power laws. This method is an efficient way to compare 
simulations with existing theories and to observe any de- 
viation from the expected scalings. As another advantage 
of compensating the energy spectrum, one notices that 
since the compensating factor increases with fcj^, it en- 
hances any possible bottleneck region occurring at large 
k^. In the numerical results that we present in this pa- 
per, no evidence of a significant bottleneck effect is found, 
which is consistent with simulations by other groups, e.g., 
[H) 13 Note that for the resolution used in our 

paper, the corresponding HD simulations would already 
produce well developed and clearly observed bottleneck 
region [36[. 

H. Locality 

The question of universality of MHD turbulence is 
closely related to the question of locality of nonlinear 
interactions. Locality loosely implies that only fluctua- 
tions of comparable length scales interact strongly with 
each other. More precisely, locality means that the prop- 
erties of the inertial interval asymptotically become inde- 
pendent of the details of the driving and the dissipation 
as the Reynolds number increases. Locality is an essen- 
tial property of hydrodynamic turbulence. In relation to 
MHD turbulence, it has been discussed in several works, 
e.g., [33-113 • Recently, Aluie and Eyink [4l| showed both 
analytically and numerically that, similarly to hydrody- 
namic turbulence, scale-locality holds in MHD turbu- 
lence. This is consistent with the fact that numerical 
simulations of strong MHD turbulence performed with 
different forcing and dissipation mechanisms produce the 
same energy spectrum, e.g., [TgI [l9l - l2ll . [23| . 

I. Integration time 

Finally, we point out that numerical simulations of 
imbalanced MHD turbulence require longer integration 
time in order to accumulate good statistics. This may 
be related to the fact that the nonlinear interaction is 
reduced in imbalanced turbulence, see subsection IIII El 
The numerical simulations presented in the next section 
indicate that for a modest imbalance of about 7 ~ 2 — 3, 
the relaxation time of the turbulence spectrum is about 
20 formally estimated large-scale dynamical times (see 
below), while about 100 dynamical times are required to 
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Run 


Resolution 


Re 










A 


512^ X 256 


2400 


4.2 X 10" 


4 


5 





B 


256^ 


900 


1.1 X 10" 


3 


10 


0.6 


C 


512^ X 256 


2200 


4.6 X 10" 


-4 


10 


0.6 


D 


1024^ X 256 


5600 


1.8 X 10" 


-4 


10 


0.6 



TABLE I: Summary of simulations of strong balanced turbu- 
lence (A) and strong imbalanced turbulence (B, C, D). 



accumulate good statistics. Numerical simulations with 
significantly shorter integration time do not produce re- 
liable results. 



A correspond to balanced turbulence, that is, Cc = 0. 
In runs labeled B, C and D, cross helicity is injected at 
the forcing scale in such a way that Gc reaches a steady 
state of ~ 0.6. We use a short time-correlated forcing, 
which is on average 1/20*'' of the Alfven time of the ex- 
cited modes, so that the energy injection rates for both 
z+ and tT only depend on the variance of the imposed 
forcing, which is controlled in our simulations. In the im- 
balanced case, field-parallel box size is optimized to reach 
the critical balance at the large scales. Except for the 
Reynolds numbers, simulations B, C, and D have the ex- 
act same parameters including the energy injection rates, 
and e~. 



IV. NUMERICAL RESULTS 

Equations ([3]) are evolved until a stationary state is 
reached, as determined by the time evolution of the to- 
tal energy of the fluctuations, (see figure [2]). A typical 
run produces over 200 snapshots; the large-scale dynamic 
time associated to the dominant large scale mode fcj^ '-^ 3 
is {L±/3)/u 

rms ^^ 
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FIG. 2: Evolution of energy for imbalanced run. The outer- 
scale dynamic time is about (Lx/3)/urms ~ 2. 

Since the background magnetic field must be strong, 
we choose So = 5 in the Vrms units, (see the discus- 
sion in [20i]). Time is normalized to the Alfven tran- 
sit time To = where L|| is the field-parallel box 
size. This time is equivalent to the perpendicular transit 
time Li_/vrms when w^ms = 1- The Reynolds number 
is defined as Re = v„ns{L^/'2,'K)/v and we have cho- 
sen the same value for the magnetic Reynolds number, 
Rm = 6rms (ij_/27r)/?7, denoting both by i? in (l3|). In 
each run, the average is performed over about 100 large- 
scale-eddy turnover times. The results are presented in 
Fig. (131). 

TableUsummarizes five representative simulations that 
incorporate all the aspects discussed in section Hill Run 
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FIG. 3: Spectra of the Elsasser fields in numerical simulations 
of strong MHD turbulence. Top two frames: balanced turbu- 
lence (run A); bottom frame: imbalanced turbulence (runs 
B(red), C(blue) and D(black)). 

The field-perpendicular energy spectrum is obtained 
by averaging the angle-integrated Fourier spectrum, 

E{k^) = 0.5(|v(ki)|2)fc^ -f 0.5(|b(ki)|2)fci, (7) 

over field-perpendicular planes in all snapshots. Figure |3l 
shows the field-perpendicular energy spectra for each run. 

Our numerical setup offers significant advantages over 
full MHD simulations, and can be already seen in sim- 
ulations of balanced turbulence, top frame in Fig. 
run A. The energy spectra approach i?*(fcj^) cx fcj^"^^^, 
in good agreement with earlier numerical findings [e.g., 
[TgI. IITI. [l9l [20j , but requiring considerably less compu- 
tational cost and producing slightly larger inertial inter- 
vals. Our most significant results are obtained for the 
imbalanced case. The bottom frame of Fig. |3l shows the 
spectra for three different Reynolds numbers (runs B, C, 
D). It is observed that the spectra are pinned at the 
dissipation scales, which supports the phenomenological 
predictions by [Ij, li^ - li^ l . We also find that the large- 
scale parts of the both spectra are practically insensitive 
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FIG. 4: Left: Cosine of the alignment angle between va and 
bA fluctuations in the guide-field perpendicular plane at scale 
A = 1/^/12 in Run A. Right: A correlated region of (counter- 
) aligned magnetic and velocity fiuctuations (red and blue vec- 
tors) at scale A — L±/12, in a plane perpendicular to the 
strong guide field, in run A. The fiuctuations are aligned pre- 
dominantly in the x direction while their directions and am- 
plitudes change predominantly in the y direction. 

to the Reynolds numbers. These two important proper- 
ties imply that as the Re numbers are further increased, 
the spectra must become progressively more parallel 
in the inertial interval. This is indeed seen in our numer- 
ical simulations (B, C, D). Our numerical simulations in- 
dicate that both spectra approach the universal scaling 
of strong MHD turbulence E^{k±) cx kj_^^^, while they 
have essentially different amplitudes and correspond to 
essentially different energy fluxes. 



V. PHENOMENOLOGICAL MODELING 

In this section we propose an explanation for the ob- 
served spectra. In the case of balanced MHD turbulence, 

3/2 

our explanation of the kj_ energy spectrum relies on 
the phenomenon of scale- dependent dynamic alignment, 
see |2ll l45l |. (General details on dynamic alignment in 
MHD turbulence can be found in, e.g., [H, li^]). Con- 
sider the eddy shown to the right of Fig. SI obtained 
from simulations. In this eddy fluctuations are aligned 
within the small angle 9x along x, while their directions 
and magnitudes change in an almost perpendicular di- 
rection, along y. In the case of strong balanced turbu- 
lence, the nonlinear interaction in such an eddy is then 
reduced by a factor 6\ for both z+ and fields, and the 
corresponding nonlinear interaction time is estimated as 
Tx l/(zj • kj.) - l/{zfk±9x). The scaling of the 
fluctuating fields is then found from the requirement of 
constant en ergy fluxes: {zf)'^/Tx = const. One can ar- 
gue [13, Hi], |45[ that the alignment angle decreases with 
scale as Ox oz A^^^, in which case the field-perpendicular 
energy spectrum is E{k±) cx kj^^^^. 

It turns out that the imbalanced simulations of previ- 
ous section can also be explained within the framework 
of the scale- dependent dynamic alignment [l^ls^]. Let 
us assume that the alignment is present in the imbal- 
anced case. Since the fields amplitudes and z~ are 
now essentially different, the alignment angles are dif- 
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FIG. 5: Sketch of dynamic alignment of magnetic and velocity 
fiuctuations in a turbulent eddy. 

ferent as well; we denote them and 6*^, see Fig. [S] 
These angles obey the important geometric constraint: 
O^z^ ~ OxZx , as is clear from Fig. [5l The depletion of 
nonlinear interaction [U IZ^ is therefore different for z+ 
and z^ fields, which makes their nonlinear interaction 
times, l/{z^k±9^), the same. The requirement of 

constant energy fluxes {z^)'^/t^ ~ = const then en- 
sures that z1^ / z^ ^ Y^e+/e^, so both fields should have 
the same scaling, although different amplitudes. We also 
note that the Reynolds numbers that take into account 
the depletion of nonlinear interaction are the same for 
both fields, i?ej ~ Xz^O^/v, and they are reduced ap- 
proximately by the factor 7 = z'^ j z~ with respect to the 
formal Reynolds number Rcx — Xvxjv- 

We note that the above relations of the kind j z'^ ~ 
y/e^~f^ hold locally for a particular domain, with given 
(positive or negative) alignment, while MHD turbulence 
consists of both positively and negatively aligned do- 
mains of various strengths. Since E'^ and e"*" are concen- 
trated mostly in positively aligned domains, while E~ 
and e~ in negatively aligned ones, the quantities {z^) 
and (e*) averaged over the global system should not a 
priori satisfy the same relations. The difference between 
the local and global quantities in MHD turbulence should 
be taken into account when one desi gns numerical tests. 

Recently, Beresnyak & Lazarian j23l [28j attempted 
to test the effects of the dynamic alignment in imbal- 
anced turbulence by measuring the relations among global 
quantities {z^) and (e^). They did not observe the local 
scaling relations and concluded that the dynamic align- 
ment proposed in [2^ [s^l is absent. According to the 
explanation given above, such a conclusion is incorrect. 
The relations between the global quantities should de- 
pend on the details of the distribution of positive and 
negative domains in the turbulent system, see the dis- 
cussion in [2^. For our present purposes we simply need 
the fact that if each domain has the same scaling of the 
fluctuating fields, the fields averaged over the whole tur- 
bulent system will have the same scaling as well. 

VI. DISCUSSION. 

We have presented a detailed numerical setup based 
on the Reduced MHD model (RMHD), and results from 
high resolution simulations of balanced and imbalanced 
MHD turbulence in steady state. We have used this nu- 
merical set up to address currently existing controversies 
regarding the spectra of imbalanced MHD turbulence. 
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The simulations are consistent with the theories and 
observatio ns p redicting same scahng for both Elsasser 
fields [H, [H^ m, [2§| and less consistent with the mod- 
els predicting different scalings for E^, [e.g., [l^]. The 
measured scaling exponent in the simulations is close to 
the —3/2 supported by phenomenological models based 
on dynamic alignment pBI . [29j . The analysis of our sim- 
ulation results may also explain somewhat puzzling nu- 
merical findings by Beresnyak and Lazarian |23l . [28j , who 
report different spectra for the Elsasser fields, and the 
intersection of the spectra rather than pinning at the dis- 
sipation scale. According to our results, the explanation 
might lie in the fact that in these simulations the imbal- 
ance was typically high, up to 7^ — (z+)^/(z^)^ ^ 1000, 
and, therefore, the universal regime of imbalanced MHD 
turbulence was not reached, see our analysis in sec- 
tion UlTll 

The phenomenology of scale-dependent dynamic align- 
ment can be applied to explain the observed spectra. In 



this phenomenology, the configuration space splits into 
eddies (domains) with highly aligned and anti-aligned 
magnetic and velocity fluctuations, where nonlinear in- 
teractions are reduced, as in left panel of Fig. 2] Even 
when the turbulence is balanced overall it still can be 
imbalanced locally, creating domains of positive and neg- 
ative cross-helicity. In each of these regions the picture 
of imbalanced turbulence presented above applies. When 
averaged over all the regions, the spectra of balanced tur- 
bulence are reproduced. 
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